#!/usr/bin/env python
import os
import numpy as np
import matplotlib.pyplot as plt

for i in range(0,1):    # Specify the number of iterations
    for j in range(0, 1):    # Specify the STRU for LAMMPS model_devi
        max_devi_f_values=[]
        for k in range(0,12): # 20 LAMMPS tasks munber for before stru
            directory = "/home/zxg/be-cu_dp_course/dpgen/run/iter.{:06d}/01.model_devi/task.{:03d}.{:06d}".format(i, j, k)
            file_path = os.path.join(directory, "model_devi.out")
            data = np.genfromtxt(file_path, skip_header=1, usecols=4)
            max_devi_f_values.append(data)
    
        # Convert the list to a numpy array
        max_devi_f_values = np.concatenate(max_devi_f_values)
    
        # Use numpy.histogram() to calculate the frequency of each value
        hist, bin_edges = np.histogram(max_devi_f_values, range=(0, 0.9), bins=100)
    
        # Normalize the frequency to percentage
        hist = hist / len(max_devi_f_values) * 100
    
        # Use matplotlib.pyplot.plot() to plot the frequency of each value
        plt.plot(bin_edges[:-1], hist, label = 'iter{:02d}__dev{:02d}'.format(i,j))
        plt.legend()
        plt.xlabel("max_devi_f eV/Å")
        plt.ylabel("Distribution %")
    
        with open(f'./iter{i:02d}_dev{j:02d}dist-max-devi-f.txt', 'w') as f:
            f.write("{:>12} {:>12}\n".format("bin_edges", "hist"))
            for x, y in zip(bin_edges[:-1], hist):
                f.write('{:>12.3f} {:>12.3f}\n'.format(x, y))
plt.savefig(f'iter{i:02d}_dev{j:02d}.png',dpi=300)
plt.clf()